Singularity formation in the Gross-Pitaevskii Equation and Collapse in BEC 
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We study a mechanism of collapse of the condensate wave function in the Gross-Pitaevskii theory 
with attractive interparticle interaction. We reformulate the Gross-Pitaevskii equation as Newton's 
equations for the particle flux and introduce a collapsing fraction of particles. We assume that 
the collapsing fraction is expelled from the condensate due to dissipation. Using this hypothesis 
we analyze the dependence of the condensate collapse on the initial conditions. We found that for 
a properly chosen negative scattering length the remnant fraction becomes larger when the initial 
aspect ratio is increased. 
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I. INTRODUCTION 

Collapse phenomena in Bose-Einstein condensates re- 
ceived a great deal of interest in the last few years 
Refs.§, §, |, |, §, g, |, §. The experimental technique 
based on a Feshbach resonance allows to tune, through 
application of an external magnetic field, the strength 
of the interatomic interaction and even switch between 
repulsive and attractive interactions Q. The realization 
of a collapse controlled by the magnetic field on a gas of 
85 Rb reveals a rich dynamical properties of the collapsing 
condensates Ref . JTo[] . 

When the number of atoms becomes sufficiently large, 
the interatomic energy overcomes the quantum pressure 
and the condensate implodes. In the course of the im- 
plosion stage the density increases in the small vicinity 
of the trap center. When the density approaches certain 
critical values a fraction of the atoms gets expelled. In a 
time period of an order of few milliseconds the conden- 
sate stabilizes. There are two observable components in 
the final stage of the condensate dynamics: remnant and 
burst atoms. The remnant atoms are the atoms which 
remain in the condensate. The burst atoms have the 
energy which is much larger then the energy of the con- 
densed atoms. There is also a fraction of atoms, which 
is not observable. This fraction is referred to as missing 
atoms. To account for dissipative losses in the collapse 
phenomena, the most part of theoretical studies reflected 
in existing literature try to modify the mean field theory 
through addition of damping terms modelling the mecha- 
nism of dissipative losses in the system. In particular, the 
results of experiment of Ref.[^0) have been analyzed by 
numerical solution of the Gross-Pitaevskii (GP) equation 
with three body recombination term Refs. JlT|, |l2| , |l3| |. A 
mechanism of losses through elastic collisions was dis- 
cussed in Rcf.Q. 

The References jll], [l2| follow the ideas of experimental 



work |To| and obtain the burst atoms fraction fitting the 
Gaussian distribution to the condensate wave function. 
The essence of the approach of Refs. [jll], [T^] is the idea 
that the final state of the condensate is comprised of the 
remnant fraction surrounded by a very dilute fraction of 
burst atoms. Hence, the only fraction which is expelled 
from the condensate is that of the missing atoms, while 
the burst atoms fraction is formed in the course of the 
condensate dynamics. 

Effective workings of a nonlinear mechanism of dissipa- 
tion in description of collapse in BEC arc inherently con- 
nected to the structure of singularity [|l5| which is formed 
in the course of the nonlinear dynamics as described by 
the Gross-Pitaevskii equation |Q. This equation reads 



2m 



V(x)^ = 0. 



(1) 



Here ^(x, t) is the wave function of the condensate, 
the external potential V^x) models the wall- less confine- 
ment (the trap), m is the mass of an individual atom, 
g = ATrh 2 m~ 1 a s is the effective interaction strength (the 
'coupling constant'), a s is the scattering length, and 
A = z2i ~§^z i s the Laplace operator. A convenient 
as well as practical choice for the confining trap is the 
paraboloidal potential V — ^ J^j w f x i 2 - Notice that 
so far very little is known about the singularity struc- 
ture of the 3D GP equation Eq.(Q) outside treatments 
of spherically symmetric case (see the monograph |l(| 
and references therein). In this work we will try to gain 
an insight into the problem through reformulation of the 
Gross-Pitaevskii equation Eq. (|l|) in a different form. Our 
approach resembles, at least formally, the dc Broglie- 
Bohm formulation of quantum mechanics [[17, 18 1. 

We introduce quantum trajectories (g-trajectories) of 
the GP equation Eq. ([!]). These trajectories satisfy a sys- 
tem of ordinary differential equations, viz. 



drj = 1 d<l> 

dt m dxi ' ' ' 



n(r),0) = rji, 



(2) 
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with i = 1,2,3 in the 3D case. Here 4> is the phase 
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defined by the representation '5 = y/pe 1 ^^ of the con- 
densate field. The fields r,(?7, i) can now be interpreted 
as trajectories of the fictitious particles with the initial 
points located at rji. These trajectories satisfy the New- 
ton equations of motion for a set of particles, 

dt 2 dxi dxi ' 
where the 'pressure potential' U(x,t) is 

U(x,t) = -^p-^A P ^+ 9P . (4) 
2m 

In addition, the Liouville formula for Eq. (Q) can be writ- 
ten in the form 

det^Jp(r,t)=po(i 7 ), (5) 

where po is the initial density of the condensate. The sys- 
tem of equations Eqs.(||), (||) is equivalent to the original 
GP equation Eq.(0). 

As was indicated above the g-trajectories are used in 
the de Broglie-Bohm causal interpretation of quantum 
mechanics Refs.jlT], pL8| . In our previous work Rcf. |hJ wc 
employed these trajectories to define singularities of the 
GP equation Eq.(jl|). Indeed, Eq.(||) immediately sug- 
gests that a set of singular points r s of the GP equa- 
tion Eq.(|l|) at which p(r s ,t*) — oo for some fixed t* in 
terms of g-trajectories, describes either the caustic or fo- 
cal points of the trajectories. Notice that in the case of 
linear Schrodinger equation the probability density p is 
always finite and caustics do not exist. 

The Newton equations Eqs.(||) define the flux of parti- 
cles. The appearance of the burst atoms in the mean field 
theory can be qualitatively understood through analysis 
of the energy distribution near the singularity in the flux 
of particles as shown in Figure [l| In this picture the peak 
with positive energy density arises due to the accelera- 
tion of the particles by the singularity. This acceleration 
is proportional to gV p and results in increase of the ki- 
netic energy of the particles surrounding the singularity. 
In Figure |l| the curves with greater peaks correspond to 
the moments of time closer to t* . The burst atoms can be 
formed only if the dissipation losses occur faster than the 
collapse domain is refilled by the flux. In this case one 
observes intermittent implosions as reported in Ref . pl| . 
The burst atoms cannot appear in the framework of the 
mean field theory if the damping parameter is chosen to 
be large enough (cf. Ref.|E|]). In this case the density 
cannot increase to sufficiently large values and the peak 
of the positive energy density is not formed. 

The dynamics of the g-trajectories which begin at the 
points located far from the origin in the ry-space can- 
not be significantly affected by possible dissipative term. 
This is because the compression time of the breathing 
mode of the harmonic oscillator, which is 1/4 of the pe- 
riod of the oscillator, is of the same order of magnitude as 
the time of the collapse including the time of dissipation. 
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FIG. 1: Energy/iV per sphere of radius \rj\ at different times 
for the spherically symmetric case. Greater peaks correspond 
to the moments of time closer to t», rj labels the initial po- 
sition of the particle in the flux of particles. Dimensionless 
units are defined in the text. 



This means that the remote trajectories cannot reach the 
collapse domain. 

In the presence of the burst atoms the condensate dy- 
namics outside the dissipation domain is a quite compli- 
cated interference of the 'outgoing' burst atoms and 'in- 
coming' atoms. Should we neglect this interference, the 
fraction of the remaining atoms could be found through 
the analysis of the particles flux falling into the singu- 
larity and described by the Newton equations Eqs.([3|). 
Notice that these equations are still valid in the presence 
of dissipative term. However in this case the relation 
Eq. (||) does not hold. 

The analysis of flux of particles falling into the singu- 
larity shows that the fraction of atoms remaining in the 
condensate is consistently constant for the wide range of 
parameters. This effect is supported by an experimen- 
tal evidence [|l0| . In this work we deal with the general 
case of asymmetric trap. In our previous work Ref. |l9f| 
we reported for the symmetric trap an effect similar to 
the described above. Our computation was based on a 
Gaussian trial wave function which was used to estimate 
a part of atoms focusing on the caustic. Here we discuss 
different types of collapse which can exist in asymmetric 
trap. 

The paper is organized as follows. In the next sec- 
tion we discuss different types of collapse which can oc- 
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cur in the asymmetric trap. In the section 3 we obtain 
the collapsing fraction and show that this fraction in the 
broad range of parameters consistently remains almost 
constant. 



II. BREATHING MODE OF THE 
CONDENSATE 

In this section we consider types of collapse which can 
be observed in an asymmetric trap. We use as the initial 
condition the Gaussian profile. For this profile 
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the harmonic oscillator equation has explicit solution, 



n(r},t) = ai{t)r]i, 



(6) 



where 



(Ti(t) = \ cos 2 (ujit) H jsm 2 (ojit), 



muji 



It seems natural, at least for Gaussian initial conditions, 
to expect that the breathing mode described by Eq. (||) 
gives the main contribution into the action 



S(v) = JdtJ d VPo (ri)£, 



where 



r 771 -2 



2mp 



(VVp) 2 



v-{, 



This provides us with the motivation to employ the 
variational principle and to seek for the solution of the 
GP equation Eq. (|l|) in the form 7^(77, t) = Ti(t)rji. Sub- 
stituting this into the action we obtain 



2 h 2 1 

Ti + UJ l Tj --r 



gN 



1 1 



(2ir)3/ 2 mi\ k T k T t 



(7) 



The density as obtained through the relation Eq. (|5|) 
reads 



p(x,t) 



N 



3 1 



2 1 2 

-X i /T i 



(8) 



Here we have included Wj in Tj and used po with Wj = 1. 

The result Eq. (^) is well-known Gaussian trial wave 
function Ref . p0[ . This function describes the breathing 
(monopole) mode of the condensate. 

In what follows we measure the time in units of 1 /u>, 
where u — (ujxUjyUJz) 1 ^ , and the distance in units of 
the the oscillator length clho = -Jft/mu). For the axial 



symmetric equation we have lo x — uj y — u r and t x — r y = 
r r . Then the system Eqs.(Q) is reduced to the following 
equations 



1 



T T Z T<! 



= 



= 



(9) 



(10) 



Here 8 = and k = 

' w a HO 

For negative k, such that \k\ is greater than certain 
critical value, the system of Eqs.(||), ( |lO| ) has singular 
solutions with r r — > as t — > . We found three different 
regimes of collapse. 

Spherically symmetric collapse is characterized by 
t z ~ TV as t ~ t* . The functions r r , t z near the sin- 
gularity are described by the power law 



T z ~ T*(t* - t) 



2/5 



^isji tj ^ 



where r* is some constant. In this regime the compression 
of the condensate is uniform in all directions. 

Radial collapse is effectively two dimensional collapse 
for which the condensate contracts only in the radial di- 
rection. The behavior of r z , r r near i* is defined by 



A. 



1/2 



(11) 
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1/4 



Here A is a constant and r r * = \/2 (y^^ 

Axial collapse is characterized by t z ~ t 2 as t ~ t*. 
In this case the compression of the condensate is faster 
in the axial direction. The function t z approaches zero 
through an oscillating regime. This regime cannot be 
described by a power law. 

The actual type of singularity depends on the ini- 
tial conditions for Eqs.(p|),([To|). To begin the analy- 
sis we specify the parameters and initial data as (3 = 
0.5325, 7,(0) = 0,f r (0) = 0,r,(0) = t z0 and r r (0) = r r0 . 
The constants t z q, JVo are a stationary solution of the sys- 
tem of Eqs. ([)[), ( p^| ) with n = n in > 0. This choice of 
parameters models experiments on Feshbach resonances 
when the scattering length is instantly changed from pos- 
itive to negative values in the initial stage of the dynam- 
ics. 

For accurate implementation of numerical simulation 
we use the following change of variables 



0M6) 



9 = In ■ 



tr.-f 



where t c is a constant. The functions v r , v z are regular in 
the interval [0, +00) if t c < t* and acquire a singularity 
at some 9 S < +00 if t c > t*. 

We also obtain A as a function of k, n in as shown 
in Figure |^. In this Figure the open domain fl\ in 
k — Ki n plane consists of two disconnected components 



FIG. 2: The 3D plot and the contour plot of A as a function 



FIG. 3: Typical behavior of u t ,v T in the axial collapse. Pa- 



rameters are |k| = 0.8, m n = 0. This figure is derived in the 
of l K l and ^' limit t c ^U-0. 



il^ 1 ' 2 -' = Q,^ (Jfi^J located on the different sides 
of the groove. The overall domain il\ corresponds to the 
radial collapse. The component il^ corresponds to the 
simple radial collapse. In this subdomain the attraction 
between atoms is strong, i.e. is large. This results in 
fast and monotonous decreasing of the Gaussian widths 
Ti, i = x,y,z. Notice that while t XiV vanish, r z remains 
finite as in Eq.(0). In the second component of £1a, i.e. 

(21 

in VL X the density p(0, t) has at least one local maximum 
in the interval (0,t»). In this regime the collapse occurs 
after a regular compression-expansion cycle. The rates 
of contraction in the radial and axial directions are alter- 
nating with each other in magnitudes. At some instances 
an intensive contraction in one direction is compensated 
by expansion in the other. Notice that the final stage 
of this process is always the contraction in the radial di- 
rection. This is why we refer to this type of collapse as 
to the radial collapse. We turn next to the description 
of the open domain complimentary to the domain fl\ in 
the k — Ki n plane. This is the groove shown in Figure 
|2|. This domain corresponds to the axial collapse. In the 
case of axial collapse the ultimate stage of the dynamics 
is always the faster contraction in the axial direction. A 
typical behavior in axial collapse is shown in Figure |[ 

The boundary separating the open domain Q\ and the 
groove, i.e. dfl\ corresponds to the spherically symmet- 
ric collapse. It is interesting to notice that even though 
the trap is asymmetric, the workings of nonlinear dynam- 
ics still allow the spherical collapse to occur. 

Intuitively, the existence of a regular cycle in dynamics 
of the breathing mode for (k, n in ) 6 suggests that 
fewer number of atoms must be removed to stabilize the 
condensate. We further discuss this idea in the next sec- 
tion where we analyze the flux of particles. 



III. THE FLUX OF PARTICLES 

The g-trajectories 7^(77, i) — Tirji derived in the previ- 
ous section focus at the symmetry axis when the time 
approaches for any r\. However, for the large |?7| this 
is not true, since in this region the behavior of the func- 
tion r(?7, t) is determined by the solution of the harmonic 
oscillator equation Eq.(||). 

The scope of analysis of equations for quantum tra- 
jectories Eqs.(||) proposed in this section attempts go 
beyond the limitations posed by the Gaussian ansatz. 
In what follows we construct an effective approxima- 
tion for the flux of particles. Our approximation effec- 
tively uses the notion of quantum trajectories. We gain 
a motivation for this approach from the analysis of the 
quasi-classical approximation for the Schrodinger equa- 
tion. To construct the quasi-classical approximation of 
the Schrodinger equation we have to solve the Newton 
equations for the trajectories of particles, i.e. to solve 
r\ = — The amplitude and the phase of the wave 
function are determined through Eqs.(||),(j|). We can 
search further corrections substituting the quasi-classical 
amplitude into the Newton equations Eqs.(|2|) with g = 0. 
Iterating this procedure and assuming that iterations 
converge we finally get the solution of the Schrodinger 
equation in the time interval [0, T] such that the func- 
tion r(?7, t) with t <G [0, T] is invertible at each step. For 
instance, in the case of harmonic oscillator and Gaussian 
initial conditions this procedure converges in the inter- 
val defined by the inequality \\ i cos(ojit) > 0. In fact 
we construct here an iteration procedure F of the form 
r (n) _ p (r'™ -1 )), where n is the number of iteration. 
These iterations start from the trajectories of the clas- 
sical harmonic oscillator = rf(rj,t) = cos^i)?^. 
Notice that the mechanism of the iteration procedure F 



is disconnected from the quasi-classical approximation. 
This procedure is equally well applicable in both linear 
and nonlinear (g ^ 0) cases. Indeed, solving Eqs.(||) in 
the nonlinear case, we can then find from Eq.(^) a next 
iteration for the function p which we again employ in 
Eqs.(||), and so forth. 

Below we introduce an effective approximation for the 
quantum trajectories as the first step of the iteration pro- 
cedure F described above. We assume that the starting 
solution Eq.(||) delivered by the variational principle is 
reasonably close to the exact solution. Using Eqs.(||) 
with the density given by equation Eq.(^J) and substitut- 
ing Eq.(@) into Eqs.(@) we obtain 



d 2 n 
dt 2 



dV. 



eff 



where 



V ef f(r,t) 
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2 ^ 



fi(0) =Vi, h(0) =0, (12) 



r? + 47r«p(r, t) /N. 
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FIG. 4: The fraction N e as a function of |«| for three different 
values of Kin- 



Notice that for the case of an axially symmetric trap it 
is sufficient to consider the behavior of the g-trajectories 
in x-z plane. This can be achieved through reduction 
r y = 0. 

At the moment of time t c close to t* we observe two 
main types of trajectories: the trajectories which go away 
from the center and the trajectories which tend to focus 
on the symmetry axis. The domain in the ?7-space corre- 
sponding to the focussing trajectories is denoted by f2 s - 
We now define the collapsing fraction of atoms by an 
integral 



podrj. 



(13) 



Here po is the initial density. It is physically plausible to 
assume that the collapsing fraction is removed from the 
condensate due to the dissipation. Then the remaining 
fraction is given by the formula 1 — N e /N. 

We analyze numerically the flux of particles falling into 
the singularity. To do this we solve Eqs.(§), ©, (jl2|). 
At the moment t = t c the focussing trajectories with rea- 
sonable accuracy get inside the domain Q c . This domain 
is restricted by the ellipsoid of revolution of an ellipse 
with the semi x and z-axes equal T z (t c ) and T r (t c ) re- 
spectively. For t c — ► t» the domain fl c reduces to the 
interval [—A, A] in the z-axis. 

We now technically define the domain fl s in the rj- 
space such that if r(?7, t c ) <G J7 C then n £ il s . The results 
of simulations are presented in Figure [|. We find that the 
fraction N e /N is consistently around 70% for large values 
of \k\ and depends neither on Ki n nor on N. More exactly 
this is true for the subdomain f2^ corresponding to the 
simple radial collapse described above. The collapsing 
fraction decreases abruptly for \k\ < |k c |, where k c is 
a certain threshold value specific to the given value of 



the parameter m n . These thresholds for three different 
values of K; n are shown in Figure 0. As was explained 

(2") 

above, below the threshold \k c \, in the domain Q\ ' and 
in the groove, we observe oscillatory dynamical patterns. 
In the course of these oscillations we observe losses in the 
collapsing fraction. For large K, n and in this region the 
fraction N e /N is found to be around 20% . 



IV. CONCLUSIONS 

In this paper we suggest a simple model Eq. (|l2] ) of the 
BEC collapse based on the dynamical properties of the 
GP equation. Our model is independent of the damp- 
ing parameter and a microscopic mechanism of the atom 
ejection. In the framework of this model it is impossible 
to correctly interpret the burst atoms fraction. However, 
the simplicity of the model allows us to study the prob- 
lem in wide range of the parameters. 

Using the notion of the flux of particles we define the 
collapsing fraction iV e as the part of atoms focusing at 
the symmetry axis. This focusing results in infinite in- 
crease of the density. Because of this we assume that 
total collapsing fraction is removed from the condensate. 
We estimate that N e /N = 0.7 for sufficiently large k 
and the remaining fraction in this case is insensitive to 
the choice of parameters. Interestingly, our results show 
that changing the initial aspect ratio by Ki n it is possible 
to move from the domain f^ 1 "* to the domain Q?' for the 
same interaction strength k. This means that the frac- 
tion of remaining atoms can be increased for the larger 
initial aspect ratio controlled by n in . In our opinion this 
effect deserves experimental investigation. 
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